#ifndef __WAVE_PARTICLES_H
#define __WAVE_PARTICLES_H

template <typename value_>
struct particle 
{
	typedef value_ value_type ;
	typedef math::vector <value_type , 2> vector2;
	typedef math::vector <value_type , 3> vector3;
	typedef vector2 point2 ;
	typedef vector3 point3 ;

	typedef particle<value_type> particle_type ;

public :
	bool is_live () const 
	{
		return pos.z >= m_min_amp ;
	}

	void set_angle (value_type slice_angle, unsigned count) 
	{
		this->a = slice_angle * count;
		this->da = slice_angle / 1.5f; 
	}
	void set_pos  (value_type x_, value_type y_, value_type amp_)
	{
		this->pos = point3(x_, y_, amp_);
	}
	void update_velocity ()
	{
		vel = point3 (m_wave_speed * cos (a), m_wave_speed * sin (a), vel.z) ;
	}
	void advent_time (value_type t)
	{
		pos += vel * t ;
	} 
	void advent_distance ( value_type d)
	{
		pos += vector3 (cos(a), sin(a), 0) * d;
	}
	// for particle spawn
	void advent_time_angle (value_type t, value_type angle)
	{
		pos += vector3 (cos(angle), sin(angle), 0) * m_wave_speed * t ;
	}
	void set_time_tag (value_type bt_tag)
	{
		bt = bt_tag ;
		st = bt + c / sin(da);
	}
	void decrease_amp ()
	{
		pos.z /= amp_dec_coeff ;
	}
	// bounding box 
	bool boundary_detect ()
	{
		return (pos.x >= 0.0f && pos.x <= 1.0f && pos.y >= 0.0f && pos.y <= 1.0f);
	}
	// reflection test
	void reflect (value_type system_time) 
	{
		if ( rx == 1 &&  tx < system_time) 
		{
			 pos.x = - pos.x;
			 a = PI -  a;
			 update_velocity();
			 rx = 0;
			 boundary_testX();
			 return ;
		}
		if ( rx == 2 &&  tx < system_time) 
		{
			 pos.x = 2.0f -  pos.x;
			 a = PI -  a;
			 update_velocity();
			 rx = 0;
			 boundary_testX();
			 return ;
		}
		if ( ry == 1 &&  ty < system_time) 
		{
			 pos.y = - pos.y;
			 a = - a;
			 update_velocity();
			 ry = 0;
			 boundary_testY();
			 return ;
		}
		if ( ry == 2 &&  ty < system_time) 
		{
			 pos.y = 2.0f -  pos.y;
			 a = - a;
			 update_velocity();
			 ry = 0;
			 boundary_testY();
			 return ;
		}	
	}
	void boundary_test ()
	{
		boundary_testX(); 
		boundary_testY();
	}

protected:
	void boundary_testX()
	{	
		value_type t =  st -  bt;
		value_type newx =  pos.x + t *  vel.x;
		if (newx < 0.0f) 
		{
			 tx =  bt - t / (newx -  pos.x) *  pos.x;
			 rx = 1;
		}
		else if (newx > 1.0f) 
		{
			 tx =  bt + t / (newx -  pos.x) * (1.0f -  pos.x);
			 rx = 2;
		}
		else  rx = 0;
	}
	void boundary_testY()
	{
		value_type t =  st -  bt;
		value_type newy =  pos.y + t *  vel.y;
		if (newy < 0.0f) 
		{
			 ty =  bt - t / (newy -  pos.y) *  pos.y;
			 ry = 1;
		}
		else if (newy > 1.0f) 
		{
			 ty =  bt + t / (newy -  pos.y) * (1.0f -  pos.y);
			 ry = 2;
		}
		else  ry = 0;
	}
private :
	bool m_live ;
public :
	static constexpr value_type get_particle_radius () {return m_particle_radius ;}
	point3 pos ;  // coordinate  (x, y), z is amplitude.  
	vector3 vel ; // velocity 

	value_type a ;    // wave angle
	value_type da, bt, st, tx, ty;

	int rx, ry;	       // reflection tag. 

private :
	static constexpr value_type m_min_amp = MIN_AMP ;
	static constexpr value_type m_wave_speed = WAVE_SPEED ;
	static constexpr value_type m_particle_radius = PARTICLE_RADIUS ;
	static constexpr value_type c = 0.25f * m_particle_radius / m_wave_speed;
	static constexpr value_type amp_dec_coeff = 2.5f ;
};


template <typename value_>
class particle_system 
{
public :
	typedef value_ value_type ;
	typedef math::vector <value_type , 2> vector2;
	typedef math::vector <value_type , 3> vector3;
	typedef vector2 point2 ;
	typedef vector3 point3 ;

	typedef particle <value_type> particle_type ;

	particle_system () : m_size (0), m_system_time(0)
	{
		this->m_particles.resize(m_max_size);
	}
	const particle_type & operator [] (unsigned i ) const {return m_particles[i] ; }
	particle_type & operator [] (unsigned i ) {return m_particles[i] ; }
	unsigned get_size () const {return m_size ;}
	unsigned get_max_size () const {return m_max_size ;}

	value_type get_system_time () const {return m_system_time ;}

	void reset () 
	{
		this->m_size = 0;
	}
	void update () 
	{
		m_system_time += m_timestep;
		recycle () ;         // collect dead particles.
		reflection_detect(); // test if particles should be bounced.
		wave_splitter() ;    // test if the particles should spawn new particles. 
	}
	void generator (value_type x, value_type y, value_type amp) 
	{
		wave_generator (x, y, amp);
	}
	
protected :
	// collect the dead particles
	void recycle ()
	{
		for (unsigned i = 0; i < m_size ; ++i)
		{
			if (false == m_particles[i].is_live())
			{
				unsigned j = m_size - 1;
				do 
				{
					if (true == m_particles[j].is_live()) break; 
					-- j; 
				}
				while (j > i) ;
				if (i != j ) 
				{
					m_particles[i] = m_particles[j] ;
				}
				m_size = j ;
			}
		}
	}
	void reflection_detect () 
	{
		for (unsigned i = 0; i < m_size; ++ i) 
			m_particles[i].reflect (m_system_time);
	}
	// generate more particles with wave propagation
	void wave_generator (value_type x, value_type y, value_type amplitude) 
	{
		value_type angle = PI360 / wave_particle_num;
		unsigned counter = m_size;
		for (unsigned i = 0; i < wave_particle_num; ++i) 
		{
			particle_type & p = m_particles[counter] ;
			p.set_pos (x, y, amplitude);
			p.set_angle (angle, i);
			p.advent_distance (wave_radius);
			if (p.boundary_detect ()) counter ++ ;else continue; 
			p.update_velocity (); 
			p.set_time_tag(m_system_time - 0.03f); // set time tag
			p.boundary_test();
		}
		m_size = counter;
	}
	void wave_splitter () 
	{
		static unsigned tick = 0; 
		if (m_size != tick)
		{
			std::cout<<"SIZE "<<m_size<<std::endl;
			std::cout<<m_particles[0].pos.z<<std::endl;
			tick = m_size;
		}
		value_type t ;
		unsigned curr_size = m_size;
		for (unsigned i = 0; i < curr_size; ++i) 
		{
			if (m_particles[i].st < m_system_time) 
			{
				particle_type curr = m_particles[i] ; 
				t = curr.st - curr.bt ;
				curr.decrease_amp();
				for (unsigned j = 0; j < 3; j++) 
				{
					unsigned k ;
					switch (j)
					{
					case 0: k = m_size ; break; 
					case 1: k = i ; break; 
					case 2: k = m_size + 1; break; 
					default : ;
					}
					particle_type & p = m_particles[k] ;
					p.set_pos (curr.pos.x, curr.pos.y, curr.pos.z); 
					if (0 == j) 
					{
						p.a = curr.a - curr.da / 2.0f ;
						p.advent_time_angle (t, curr.a - curr.da / 1.5f);
					}
					else if (2 == j) 
					{
						p.a = curr.a + curr.da / 2.0f ;
						p.advent_time_angle (t, curr.a + curr.da / 1.5f);
					}
					else p.advent_time (t) ; 

					p.update_velocity(); 
					p.da = curr.da / 3.0f;
					p.set_time_tag(curr.st); //st->bt
					p.boundary_test();
				}
				m_size += 2;
			}
		}		
	}


private:
	unsigned m_size ;
	value_type m_system_time ;

	std::vector <particle_type > m_particles; 
	static constexpr unsigned m_max_size = MAX_PARTICLE_NUMBER ;
	static constexpr value_type m_timestep = TIME_STEP ;
	static constexpr unsigned wave_particle_num = 10; 
	static constexpr value_type wave_radius = particle_type::get_particle_radius() * 0.05f / sin(4.5f * PI / wave_particle_num);
};


#endif 
